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Abstract 



h^ ■ Blood vessel networks form by spontaneous aggregation of individual 

^H , cells migrating toward vascularization sites (vasculogenesis). A success- 

(~| ■ ful theoretical model of two dimensional experimental vasculogenesis has 

been recently proposed, showing the relevance of percolation concepts and 
of cell cross-talk (chemotactic autocrine loop) to the understanding of this 
self-aggregation process. Here we study the natural 3D extension of the 
computational model proposed earlier, which is relevant for the investiga- 
^v^ . tion of the genuinely threedimensional process of vasculogenesis in verte- 

^ ' brate embryos. The computational model is based on a multidimensional 

V^Q , Burgers equation coupled with a reaction diffusion equation for a chemo- 

c7^ ' tactic factor and a mass conservation law. The numerical approximation 

^^D ' of the computational model is obtained by high order relaxed schemes. 

\l , Space and time discretization are performed by using TVD schemes and, 

respectively, IMEX schemes. Due to the computational costs of realistic 
simulations, we have implemented the numerical algorithm on a cluster 
for parallel computation. Starting from initial conditions mimicking the 
f^ ' experimentally observed ones, numerical simulations produce network-like 

"ti I structures qualitatively similar to those observed in the early stages of in 

a vivo vasculogenesis. We develop the computation of critical percolative 

indices as a robust measure of the network geometry as a first step towards 
^ , the comparison of computational and experimental data. 



b ! 1 Introduction 



In recent years, biologists have collected many qualitative and quantitative data 
on the behavior of microscopic components of living beings. We are, however, 
still far from understanding in detail how these microscopic components interact 
to build functions which are essential for life. A problem of particular interest 
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which has been extensively investigated is the formation of patterns in biological 
tissues J2| . Such patterns often show self-similarity and scaling laws JlBj similar 
to those emerging in the physics of phase transitions ;26|. 

The vascular network 01 ^^ is a typical example of natural structure char- 
acterized by non trivial scaling laws. In recent years many experimental in- 
vestigations have been performed on the mechanism of blood vessel formation 
both in living beings and in in vitro experiments. Vascular networks form 
by spontaneous aggregation of individual cells travelling toward vascularization 
sites (vasculogenesis) . A successful theoretical model of two dimensional ex- 
perimental vasculogenesis has been recently proposed, showing the relevance of 
percolation concepts and of cell cross-talk (chemotactic autocrine loop) to the 
understanding of this self-aggregation process. 

Theoretical and computational modelling is useful in testing biological hy- 
potheses in order to explain which kind of coordinated dynamics gives origin 
to the observed highly structured tissue patterns. One can develop computa- 
tional models based on simple dynamical principles and test whether they are 
able to reproduce the experimentally observed features. If the basic dynamical 
principles are correctly chosen, computational experiments allow to observe the 
emergence of complex structures from a multiplicity of interactions following 
simple rules. 

Apart from the purely theoretical interest, reproducing biological dynamics 
by computational models allows to identify those biochemical and biophysical 
parameters which are the most important in driving the process. This way, 
computational models can produce a deeper understanding of biological mech- 
anisms, which in principle may end up having relevant practical consequences. 
It is worth noticing here that a complete understanding of the vascularization 
process is possible only if it is considered in its natural threedimensional setting 

(mia). 

In this paper we illustrate computational results regarding the simulation of 
vascular network formation in a threedimensional environment. We consider the 
threedimensional version of the model proposed in JOIEBI- The model is based 
on a Burgers-like equation, a well studied paradigm in the theory of pattern 
formation, integrated with a feedback term describing the chemotactic autocrine 
loop. The numerical evolution of the computational model starting from initial 
conditions mimicking the experimentally observed ones produces network-like 
structures qualitatively similar to those observed in the early stages of in vivo 
vasculogenesis. 

Since in the long run we are interested in developing quantitative comparison 
between experimental data and theoretical model, we start by selecting a set of 
observable quantities providing robust quantitative information on the network 
geometry. The lesson learned from the study of twodimensional vasculogenesis is 
that percolative exponents jJZ] are an interesting set of such observables, so we 
test the computation of percolative exponents on simulated network structures. 

A thorough quantitative comparison of the geometrical properties of experi- 
mental and computational network structures will become possible as soon as an 
adequate amount of experimental data, allowing proper statistical computation. 



will become available. 

The paper is organized as follows. Section 2 summarizes some background 
knowledge on the biological problem of vascular network formation. Section 
3 is a short review of the properties of the model introduced in JOIEBI- In 
Section 4 the numerical approximation technique for the model is described. In 
Section 5 we describe the qualitative properties of simulated network structures 
and present the results of the computation of the exponents of the percolative 
transition. Finally, in the Conclusions, we point out at predictable developments 
of our research. 



2 Biological background 

To supply tissues with nutrients in an optimal way, vertebrates have devel- 
oped a hierarchical vascular system which terminates in a network of size- 
invariant units, i.e. capillaries. Capillary networks characterized by intercap- 
illary distances ranging from 50 to 300 /im are essential for optimal metabolic 
exchange [llj . 

Capillaries are made of endothelial cells. Their growth is essentially driven 
by two processes: vasculogenesis and angiogenesis 0. Vasculogenesis consists of 
local differentiation of precursor cells to endothehal ones, that assemble into a 
vascular network by directed migration and cohesion. Angiogenesis is essentially 
characterized by sprouting of novel structures and their remodelling. 

In twodimensional assays, the process of formation of a vascular network 
starting from randomly seeded cells can be accurately tracked by videomi- 
croscopy ^^ and it is observed to proceed along three main stages: i) migration 
and early network formation, ii) network remodelling and iii) differentiation in 
tubular structures. During the first phase, which is the most important for de- 
termining the final geometrical properties of the structures, cells migrate over 
distances which are an order of magnitude larger than their radius and aggre- 
gate when they adhere with one of their neighbours. An accurate statistics of 
individual cells trajectories has been presented in ^^, showing that, in the first 
stage of the dynamics, cell motion has marked directional persistence, pointing 
toward zones of higher cell concentration. This indicates that cells communicate 
through the emission of soluble chemical factors that diffuse (and degrade) in the 
surrounding medium, moving toward the gradients of this chemical field. Cells 
behave Hke not-directly interacting particles, the interaction being mediated by 
the release of soluble chemotactic factors. Their dynamics is well reproduced 
by the theoretical model proposed in jlOj . 

The lessons learned from the study of in vitro vasculogenesis is thus that 
the formation of experimentally observed structures can be explained as the 
consequence of cell motility and of cell cross-talk mediated by the exchange of 
soluble chemical factors (chemotactic autocrine loop). The theoretical model 
also shows that the main factors determining the qualitative properties of the 
observed vascular structures are the available cell density and the diffusivity and 
half-life of the soluble chemical exchanged. It seems that only the dynamical 



rules followed by the individual cell are actually encoded in the genes. The 
interplay of these simple dynamical rules with the geometrical and physical 
properties of the environment produces the highly structured final result. 

At the moment, no direct observation of the chemotactic autocrine loop 
regulating vascular network formation is available, although several indirect 
biochemical observations point to it, so, the main evidence in this sense still 
comes from the theoretical analysis of computational models. 

Several major developments in threedimensional cell culture and in cell and 
tissue imaging allow today to observe in real time the mechanisms of cell mi- 
gration and aggregation in threedimensional settings |9ll2T]. 

In the embryo, endothelial cells are produced and migrate in a threedi- 
mensional scaffold, the extracellular matrix. Migration is actually performed 
through a series of biochemical processes, such as sensing of chemotactic gra- 
dients, and of mechanical operations, such as extensions, contractions, and de- 
grading of the extracellular matrix along the way. 

The evidence provided by twodimensional experimental vasculogenesis sug- 
gests that cell motion can be directed by an autocrine loop of soluble chemoat- 
tractant factors also in the real threedimensional environment. 

As a sample of typical vascular structures that are observe in a threedimen- 
sional setting in the early stages of development of a living being, we include 
here (750 /im)^ images of chick embryo brain at different development stages 
(Fig.QJ. At an early stage (about 52-64 hours) one observes a typical immature 
vascular network formed by vasculogenesis and characterized by a high density 
of similar blood vessels (Fig. H^). At the next stage (70-72 hours) we observe 
initial remodelling of the vascular network (Figs. nj3,C). RemodeHng becomes 
more evident when the embryo is 5 days old, when blood vessels are organized 
in a mature, hierarchically organized vascular tree (Fig.^D). 

3 Mathematical model of blood vessel growth 

The multidimensional Burgers' equation is a well-known paradigm in the study 
of pattern formation. It gives a coarse grained hydrodynamic description of 
the motion of independent agents performing rectilinear motion and interacting 
only at very short ranges. These equations have been utilized to describe the 
emergence of structured patterns in many different physical settings (see e.g. 
OEI)- In the early stages of dynamics, each particle moves with a constant 
velocity, given by a random statistical distribution. This motion gives rise to 
intersection of trajectories and formation of shock waves. After the birth of these 
local singularities regions of high density grow and form a peculiar network-like 
structure. The main feature of this structure is the existence of comparatively 
thin layers and filaments of high density that separate large low-density regions. 
In order to study and identify the factors influencing blood vessel forma- 
tion one has to take into account evidence suggesting that cells do not behave 
as independent agents, but rather exchange information in the form of soluble 
chemical factors. This leads to the model proposed by Gamba et al. in JHI 




Figure 1: Vascular networks formed by vasculogenesis in chick embryo brain, at 
various stages of development, classified according to Hamilton and Hamburger 
(HH). A: HH stage 17, corresponding to 52-64 hours; B: HH stage 20 (70-72 
hours); C,D: HH stage 26 (5 days). 



and Serini et al. in [2SI- The model describes the motion of a fluid of ran- 
domly seeded independent particles which communicate through emission and 
absorption of a soluble factor and move toward its concentration gradients. 

3.1 Model equations 

The cell population is described by a continuous density n(x, i), where x £ R*^ 
{d = 2,3) is the space variable, and t > is the time variable. The population 
density moves with velocities v(x, t), that are stimulated by chemical gradients 
of a soluble factor. The chemoattractant soluble factor is described by a scalar 
chemical concentration field c(x, i). It is supposed to be released by the cells, 
diffuse, and degrade in a finite time, in agreement with experimental observa- 
tions. 

The dynamics of the cell density can be described by coupling three equa- 
tions. The first one is the mass conservation law for cell matter, which expresses 
the conservation of the number of cells. The second one is a momentum bal- 
ance law that takes into account the phenomenological chemotactic force, the 
dissipation by interaction with the substrate, the phenomenon of cell directional 
persistency along their trajectories and a term implementing an excluded vol- 
ume constraint jlOl E]. Finally there is a reaction-diffusion equation for the 
production, degradation and diffusion of the concentration of the chemotactic 
factor. One then has the following system: 



+ V • (nv) = (la) 

-H V Vv = ^(c)Vc - V0(n) - /3(c)v (lb) 



dn 
'dt 
dv 

'm 

dc ^ . , , c _ , 

— ^ DAc + aic)n (Ic) 

Ot T 

where /i measures the cell response to the chemotactic factor, while D and r are 
respectively the diffusion coefficient and the characteristic degradation time of 
the soluble chemoattractant. The function a determines the rate of release of 
the chemical factor. The friction term — /?v mimics the dissipative interaction 
of the cells with the extracellular matrix. 

A simple model can be obtained by assuming that the cell sensitivity /x, the 
rate of release of the chemoattractant a and the friction coefficient f3 are con- 
stant. A more realistic description may be obtained including saturation effects 
as functional dependencies of the aforementioned coefficients on the concentra- 
tion c. 

The term V0(n) is a density dependent pressure term, where (j){n) is zero for 
low densities, and increases for densities above a suitable threshold. This pres- 
sure is a phenomenological term which models short range interaction between 
cells and the fact that cells do not interpenetrate. 

We observe that, at low density n and for small chemoattractive gradients, 
ijlbb l is an inviscid Burgers' equation for the velocity field v j^, coupled to the 
standard reaction-diffusion equation ifTch l and the mass conservaton law iflah ). 



Since in the early stages of development almost all intraembryonic mesoder- 
mal tissues contain migrating endothelial precursors, we use initial conditions 
representing a randomly scattered distribution of cells, i.e., we throw an assigned 
number of cells in random positions inside the cubic box, with zero initial veloc- 
ities and zero initial concentration of the soluble factor, with a single cell given 
initially by a Gaussian bump of width a of the order of the average cell radius 
(~ 15/im) and unitary weight in the integrated cell density field n. 

In order to model the fact that closely packed cells resist to compression, 
a phenomenological, density dependent, pressure V(j>{n) acting only when cells 
become close enough to each other is introduced. The potential (j> has to be 
monotonically increasing and constant for n < uq where no is the close-packing 
density. Our simulations suggest that the exact functional form of 0(n) is not 
relevant. For simplicity we choose 

I n < no 

3.2 Parameter values 

Fourier analysis of Eq. ifTcIl with constant parameters and in the fast diffusion 
approximation dc/dt = suggests that starting from the aformentioned initial 
conditions, equation iQJ should develop network patterns characterized by a 
typical length scale ro = VDt, which is the effective range of the interaction 
mediated by soluble factors. As a matter of fact, Fourier components Ck of the 
chemical field are related to the Fourier components of the density field fik by 

the relation 

arfik 

^"^ Drk^ + r 

This means that in equation (QJ wavelengths of the field n of order rg are 
amplified, while wavelengths A ^ rg or A <C tq are suppressed. 

Initial conditions introduce in the problem a typical length scale given by the 
average cell-cell distance L/^/N, where L is the system size and N the particle 
number. The dynamics, filtering wavelengths |HI, rearranges the matter and 
forms a network characterized by the typical length scale tq. 

It is interesting to check the compatibility of the theoretical prediction with 
physical data. From available experimental results P^ it is known that the order 
of magnitude of the diffusion coefficient for major angiogenic growth factors is 
D — 10~^ cm^ s~^. In the experimental conditions that were considered in JHI 
the half life of soluble factors is 64 ± 7 min. This gives tq ^ 200 ^m, a value in 
good agreement with experimental observations. 

3.3 Lower dimensional models 

In order to get some intuition about the typical system dynamics, we exploit the 
ID version of model ^ to simulate the "collision" of two cells. For small values 
of Bp and sufficiently high Cpin ^, the two bumps merge into a single one (see 



Fig. 121 left) which appears to be stationary, as suggested also by the graphs of 
the kinetic energy and of the momentum of inertia (Fig. El top). On the other 
hand a less smooth onset of pressure obtained with larger Bp or smaller Cp 
leads to forces overcoming the chemical attractive ones, making the two bumps 
bounce back (Fig. |2lright. Fig. Elbottom). We observe that the better dynamics 
from the biological point of view is the first behavior with two bump coalescing. 

Biological observations suggest that the dynamics of cell changes when they 
estabHsh cell-cell contacts. It is reasonable to suppose that a different genetic 
program is activated at this moment, disabhng cell motility. We therefore switch 
off cell motility as soon as the cell concentration, signalled by chemoattractant 
emission, reaches a given threshold. In this way the computational system is 
guaranteed to reach a stationary state. 

These effects can be taken into account using a non-constant sensitivity ^(c), 
a non-Hnear emission rate a(c), or a variable friction coefficient (3{c). We choose 
a threshold cq and functions of the form 

fi{c) = /^o[l — tanh(c — cq)] (3a) 

a{c) = ao[l ~ tanh(c — cq)] (3b) 

I3{c) = l3o[l + tanh(c - cq)] (3c) 

The effect of the first two terms is that the sensitivity of the cells and their 
chemoattractant production is strongly damped when the concentration c reaches 
the threshold cq . We did not observe a significant dependence on the exact form 
of the damping function, provided that it approximates a step function that is 
nonzero only when c < cq. 

/3(c), on the other hand has the effect of turning on a strong friction term 
at locations of high chemoattractant concentration. We performed several tests 
and observed that the different choices are approximately equivalent in freezing 
the system into a network-like stationary state. 



4 Numerical methods 

Our scheme is based on a suitable relaxation approximation ^j of the mass 
conservation law (|la|l and the multidimensional Burgers equation lllb|l coupled 
with a second order finite-differences method for the reaction-diffusion equation 
ifTcIl of the chemotactic factor. We point out that also for the last equation ltTc|l 
we could consider a relaxation approximation X3,s ;1,8| in order to deal with the 
system (QJ in an uniform way, but we prefer to adopt here a simpler approach. 
We first briefiy review an extension of the approach proposed by Jin and Xin 
in 231 for a scalar conservation law to the case when a source term is present 

I + !;/(.) ^9(»). (4) 

Introducing an auxiliary variable j that plays the role of a physical fiux we 
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Figure 2: Bump coalescence driven by chemotactic force and pressure. In the 
first three rows the density and velocity fields at subsequent instants of time are 
shown. In the last row we show the time evolution of the kinetic energy and 
of the momentum of inertia. Left column: Cp — 3 and Bp — 10^^, leading to 
bump coalescence. Right column: Cp =§2 and Bp — 10^^, leading to undesired 
rebound of the two bumps. 
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Figure 3: Time evolution of the kinetic energy and of the momentum of inertia. 
Top: Cp = 3 and Bp = 10~^, leading to bump coalescence. Bottom: Cp = 2 
and Bp — 10"^, leading to undesired rebound of the two bumps. 
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consider the following relaxation system: 



du dj , , 
dj du 1 



f{u)l 



(5a) 
(5b) 



where e is a small positive parameter, called relaxation time, and a is a suitable 
positive constant. Formally, Chapman-Enskog expansion justifies the agreement 
of the solutions of the relaxation system with the solutions of the equation 



du 
'di 



d_ 

dx 



f{u) = g{u) 



d_ 
dx 



{a-f'iuf 



du 
dx 



(6) 



which is a first order approximation of the original balance law I^J • 

It is also clear that © is dissipative, provided that the subcharacteristic 
condition a > f'{u)^ is satisfied. We would expect that appropriate numerical 
discretization of the relaxation system Q yields accurate approximation to the 
original equation Q when the relaxation parameter e is sufficiently small. 

In view of its numerical approximation, the main advantage of the relaxation 
system Q over the original equation 1^ lies in the linear structure of the char- 
acteristic fields and in the localized low order term and this avoids the use of 
time consuming Riemann solvers. Moreover, proper implicit time discretization 
can be exploited to overcome the stability constraints due to the stiffness and 
to avoid the use of non-linear solvers. 

We observe that system ^ is in the form 



-^ +diy fiz)^g{z) + -h{z) 
dt e 



(7) 



where z = {u,j)'^, f{z) = {j,au)'^, g{z) = {g{u),0)'^ and h{z) = (0,j - f{u))'^. 
When £ is small, the presence of both non-stiff and stiff terms, suggests the use 
of IMEX schemes ji [H EO] . 

Assume for simpHcity to adopt a uniform time step At and denote with z" 
the numerical approximation at time i„ = nAt, for n = 0, 1, ... In our case a 
i^-stages IMEX scheme reads 



z"+i = z" 



At'^k 



i=l 



g(z«)+,(z«) 



^^M,(z«) 



where the stage values are computed as 



^) = z"-At^a,,fe 



fc=i 



^•^'■-.(fe)^ , „rr(fe)^ 



dx 



(zW)-h5(^^'-0 



At 



y^a».fc/t(^ 



(k)) 



fc=l 



Here {aik,hi) and {aik,bi) are a pair of Butcher's tableaux of, respectively, a 
diagonally implicit and an explicit Runge-Kutta schemes. 
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In this work we use the so-called relaxed schemes, that are obtained letting 
£ ^ in the numerical scheme for ^. For these the first stage 



,(1) 
7(1) 



u 



H aiih 

e 



,•(1) 



becomes 



^(1) = u" 



then it reduces to h{z^'^^) = 0. While the second stage, i = 2, reads 



.(2) - ,« 



z" - Ata2,i 



g(z«)+,(zW) 



=0 



which implies that h{z^'^>) — 0. 

Summarizing, the relaxed scheme yields an alternation of relaxation steps 



h{z^'^) = 



I.e. J 



ii) - 



fiu^'^) 



and transport steps where we advance for time cii^k^t 

dz 



dt 



+ div/(z) = g{z) 



with initial data z = z^^^ retain only the first component and assign it to u^^^^\ 

Finally the value of m"'+^ is computed as m" + ^ biU^'^\ 

In order to obtain a relaxation approximation of the first and second equation 
of (QJ we rewrite them in conservative form, introducing the moment p(x, f) = 
n(x, t)v(x, f): 



dn 
'dt 
dp 
dt 



+ V-p = 

+ V • (nv (E)\-) — TifiWc — nW(j>{n) — (3p 



(8a) 
(8b) 



Introducing the variable u — {n, p)-^ and the auxiliary fiux w, the relaxation 
system reads 



du ^ ^, , 

— - + V • w = G(u,w,c) 
dt 



dt 



1. 

e 



^(u)) 



(9a) 
(9b) 



where G(u, Vir, c) — {0,nfi^c — n\I(f){n) ~ /3p)"^, F{u.) ~ (p,nv ® v) and A 
is a suitable diagonal matrix whose positive diagonal elements verify a sub- 
characteristic condition. As we previously remarked, our relaxed scheme takes 
alternatively an implicit step and an explicit one: the explicit step involves the 
computation of the fiux V • w and the evaluation of the non stiff source term G. 
In particular we compute Vc and V0(n) using a second order difference scheme. 
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duj 
dt 


+ 


I (w, + l/2 - 


w,- 


dt 


+ 


^ 


(u, + l/2 


Uj 



In the following we describe for simplicity the fully discrete scheme in one 
dimensional case. We introduce the spatial grid points Xj with uniform mesh 
width h — Xj+i —Xj. As usual, we denote by u" the approximate cell average of 
a quantity u in the cell [xj-i/2,Xj^i/2] at time i„ and by m?+i/2 the approximate 
point value of u at x = 0:^+1/2 and t — tn- A spatial discretization to (EJ in 
conservation form can be written as 

1/2) =G(uj-,Wj,Cj) (10a) 

_i/2)--i(w,-f(u,)). (10b) 

In order to compute the numerical fluxes Wj±i/2, we consider the characteris- 
tic variables w ± A^''^\i that travel with constant velocities ±^^'^, and so the 
semidiscrete system becomes diagonal. Now we have to apply a numerical ap- 
proximation to w ± A^'^u. A first idea is to apply a ENO or WENO approach 
(see e.g. ^j), to build an high order reconstruction, coupled with a suitable 
IMEX scheme. The drawback is the high computational costs, especially in 
a multidimensional framework. Therefore we chose a suitable compromise be- 
tween the computational cost and the accuracy, using a second order TVD 
scheme. The numerical fiux that we use is obtained coupling an upwind scheme 
and the Lax-Wendroff method by a non linear flux limiter J^I- Namely the high 
order flux F{U) for a generic variable U consists of the low order term Fl[U) 
plus a second order correction Fh[U): 

F{U) = Fl{U) + ^iU)iFHiU) - Fl{U)) 

where ^ is the flux limiter. When the data U is smooth, then ^(C/) should be 
near 1, while near a discontinuity we want ^{U) close to 0. The idea consists 
in the selection of a high order flux Fh that works well in smooth regions and 
of a low order flux Fl which behaves well near discontinuities. 

In our schemes we considered the upwind scheme as a low order flux for the 
characteristic variables 

Fl{{^+A^'^vl),+^/2) = (w+Ai/2^)„ Fl((w-A1/2u),+i/2) - (w-Ai/2^),+i 

and the Lax-Wendroff scheme as a high order flux for the same variables 

Fff ((w ± Ai/2u)^.+i/2) - 4^((w ± A^/^n),+, + (w ± A^/'^n),) 



\A 



1/2 



2 



((w ± A^''^vl)j+^ - (w ± Ai/2u)j-) 



where A = /S.t/h (we advance of one time step). 
Letting 

^ / (W ± Al/2u)^" - (W ± yll/2„)n_^ ' 



±1 







^ Wv^f± Al/2u)" ^ - (W± Al/2u)« 



the fully discrete scheme for the variable u using Euler method to advance in 
time is the following 

Uj — Uj H 2~vUj + l ~ ^Uj + ^j-1) ~ 2\^3 + l ~ ^3-11 

A.I^AAl^/ + , + , - _ -^ 
^I 4 >^ ^3 ^ '^i-i + '^i+l *J ''' 
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with 



f = i(±Ai/2u7±i + w^'ii T A'/'u- - w«)v|/(e±). 



After the substitution of the relaxing step we get 



(11) 



,n+l 



U 

At 



n I XA^ 






(-4 + 4-1 + ^7+1 -.^7), 



where s* is obtained from ljll|l letting w = -F(u). The scheme can be put 
in a conservative form and it is possible to prove its consistency by standard 
technique ,173. In order to prove a TVD stability, we write 






where 



D" 



(1 - c] - d;o(u5'+i - up + C7_i(u7 - u»_i) 

+D7+i(u^^_2 - U^+i) + Ej+i/2, 



2^ + "7+1-"" 

A /^l/2 _ f^K+l)-FK) 

2 V ";+!-"" 



(12) 



^;Vi/2 = Ati-^ 



o'^i+2 2s^-_|_i 



4-^;+i+24-^;^ij 



where we notice that C and D are non negative. 

The coefHcient E can be written in terms of C and D, in fact 



si = 



-cr*(e+)(u;vi " u7)), ^7 = - ^-^_^^(&-){vq - u-_j). 



XAx * ^ » " J+^ J"' J xAx 

We can rewrite II12|I in the following form 



,n+l 
ij + 1 



n+1 



(u7+i - u7) [(1 - C] - D7) + (1 - AAi/2)(D7v|/7^i + C7*+)] 



C" 



l-A^i/^ 



+K 



i-1 
D 



j+i 



i^(Dpi*7 + cpi*jii: 

i-^(D;^iVi/^2 + CPi*^i; 



(13) 
It's easy to see that under the CFL condition ||A-\/max{ai}| < 1, where Oi are 
the positive diagonal elements of the matrix A, and using the fact that the flux 

limiter verifles 

*(■©) 
< — ^ < 2, < *(0) < 2, 







we have 



(i-C7-Dp + (i-AAi/2)(D7*7+i + C7*+) > 
C7_i-i^Y^(D7_i*7 + C7_iVi/+_i) > 

D7fi-^-^(D."+iVl/7+2 + C" *+ ) > 



and so we can deduce that our scheme is TVD stable from Harten's Theorem 

m 
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Scaling on the cluster ULISSE 
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Figure 4: Scaling of the 3D algorithm on the ULISSE cluster. Dots repre- 
sent execution time (s) and asterisks the number of Mflops/s for our numerical 
algorithm. Dashed and dash-dot lines are linear interpolations 

In the case of multidimensions, a similar discretization can be applied to each 
space dimension |14[ll3ini?] . Then, since the structure of the multidimensional 
relaxation system is similar to the ID system, the numerical implementation 
for higher dimensional problems, based on additive dimensional splitting, is not 
much harder than for ID problems. 

For our threedimensional problem the computational cost is quite high and 
can be reduced using parallel computing: the semilinearity of relaxation systems, 
together with our suitably chosen discretizations, provides parallel algorithm 
with almost optimal scaling properties. In particular the domain is divided 
in smaller subdomains and each subdomain is assigned to a processor. The 
computations of all non linear terms involve only pointwise evaluations and it 
is easy to perform these tasks in a local way. Only point near the interfaces 
between different subdomain need to be communicated in the transport step. 
We implemented these algorithm on a high performance cluster for parallel 
computation installed at the Department of Mathematics of the University of 
Milano (http://cluster.mat.unimi.it/). The scaling properties of the algorithm 
are shown in Fig. Eland are essentially due to the exclusive use of matrix- vectors 
operations and to the avoidance of solvers for linear or non-linear systems. 

5 Numerical results 

We perform threedimensional numerical simulations of model (QJ on a cubic box 
with side of length L = 1mm, with periodic boundary conditions. The initial 
condition is assigned in the form of a set of gaussian bumps with a = 15/im 
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1 



Figure 5: Initial state of a numerical simulation with 2500cells/mm . The 
colorbar on the right is referred to the coloring of the cross sections. The red 
three-dimensional isosurface corresponds to the black contour Hues in the cross 
sections 



scattered in the cube with uniform probability and having zero initial velocity. 

Biochemical data ^^l suggest the values D = 10~'^mm^/s and r = 4000 s 
for the diffusion constant and the chemoattractant decay rate. We fix the other 
constant parameters by dimensional analysis and fitting to the characteristic 
scales of the biological system. In particular, we choose: hq — 10~^^mm^/s^, 
a = ls~^, /? = 10~'^s~^. For the coefficients in the expression (j^J of the pressure 
function (j) we take no = 1.0, Cp = 3 and Bp = lO^-^. 

Very fine grids have to be used in order to resolve the details of the n{x,t) 
field, which may contain hundreds of small bumps, each representing a single 
cell. Since each cell has radius a = 15/im, one needs a grid spacing such that 
Ax < 10/im and therefore grids of at least 100^ cells for a cubic domain of 1mm 
side. 
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Figure 6: Transient state of the evolution of the initial state depicted in Fig. 
according to model Q . The initial formation of network- like structures is 
observed 
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Figure 7: Stationary state of the evolution of the states depicted in Figs. 
and El according to model (QJ . Well developed threedimensional network-like 
structures are observed 
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We performed numerical simulations with varying initial average cell density 
fi. We observed that the initially randomly distributed cells coalesce forming 
elongated structures and evolve towards a stationary state mimicking the ge- 
ometry of a blood vessel network in the early stages of formation. 

We assigned fi in the range 2100 — 3500 cells/mm"^ and performed 10 to 15 
runs for each density value with a 128 x 128 x 128 grid on a biological system 
of 1 mm^. The characteristic lengths and geometric properties of the stationary 
state depend on n and we observed a percolative phase transition similar to the 
one described in pHI for the twodimensional case. 

5.1 Analysis of the percolative phase transition 

In experimental blood vessel formation it has been shown that a percolative 
transition is observed, by varying the initial cell density. For low cell densi- 
ties only isolated clusters of endothelial cells are observed, while for very high 
densities cells fill the whole available space. In between these two extreme be- 
haviours, close to a critical cell density Uc, one observes the formation of critical 
percolating clusters connecting opposite sides of the domain, characterized by 
well defined scaling laws and exponents. These exponents are known not to 
depend on the microscopic details of the process while their values characterize 
different classes of aggregation dynamics. 

The purely geometric problem of percolation is actually one of the simplest 
phase transitions occurring in nature. Many percolative models show a second 
order phase transition in correspondence to a critical value ric, i-e. the proba- 
bility n of observing an infinite, percolating cluster is for n < ric and 1 for 
n > ric [23 • The phase transition can be studied by focusing on the values of an 
order parameter, i.e. an observable quantity that is zero before the transition 
and takes on values of order 1 after it. In a percolation problem the natural 
order parameter is the probability P that a randomly chosen site belongs to the 
infinite cluster (on finite grids, the infinite cluster is substituted by the largest 
one). 

In the vicinity of the critical density ric the geometric properties of clusters 
show a pecuHar scaHng behavior. For instance, in a system of Hnear finite size 
L, the probability of percolation n(n, L), defined empirically as the fraction 
of computational experiments that produce a percolating cluster, is actually a 
function of the combination {n — Uc) L^'^ , where i^ is a universal exponent. 

In a neighborhood of the critical point and on a system of finite size L, the 
following finite size scaHng relations are also observed: 

n(n,L) - n[(n-nc)L^/''] (14) 

There are two main reasons to study percolation in relation to vascular net- 
work formation: ( i) percolation is a fundamental property for vascular networks: 
blood should have the possibility to travel through the whole vascular network 
to carry nutrients to tissues; (m) critical exponents are robust observables char- 
acterizing the aggregation dynamics. 
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A rather complete characterization of percolative exponents in the two- 
dimensional case has been provided in |10j . 

As a first step in the study of the more realistic threedimensional case, we 
compute the exponent v characterizing the structures produced by the model 
dynamics (QJ with varying initial cell density. 

To this aim, extensive numerical simulation of system (QJ were performed 
using lattice sizes L = 1,0.78,0.62,0.5 mm, with different values of the initial 
density n. For each point 10 to 15 reaHzations of the system of size 1mm were 
computed, depending on the proximity to the critical point. 

The continuous density at final time n(x) was then mapped to a set of 
occupied and empty sites by choosing a threshold tlq. Each region of adjacent 
occupied sites (cluster) was marked with a different index. The percolation 
probability 11 for each set of realizations was then measured. In Fig. |S| we 
show clusters obtained in a box with L = 0.5 mm with n = 3100. The largest 
percolating cluster is shown in red, together with some other smaller clusters 
shown in different colors. 

Using relation (I14II . we estimate the position of the critical point Uc and 
the value of the critical exponent ly. The data for different box side length 
and initial density should lie on a single curve after rescaling the densities as 
n = {n — nc)L^^'^. For fixed Uc and v we rescale n and fit the data with 
a logistic curve, then compute the distance of the data from the curve. The 
squared distance is minimized to obtain estimates for ric and u. 

Using no — 0.35 we obtain ric = 2658 and v = 0.84. This latter value is com- 
patible with the known value 0.88 for random percolation in three dimensions 

EH. 



6 Conclusions 

We have exposed results on the numerical simulation of vascular network for- 
mation in a threedimensional setting. 

We have used the threedimensional version of the equations proposed in 
jlOl 12^1 as a computational model. Evolution starting from initial conditions 
mimicking the experimentally observed ones produce network-like structures 
qualitatively similar to those observed in the early stages of in vivo vasculoge- 
nesis. 

As a starting point towards a quantitative comparison between experimental 
data and the theoretical model we nedd to select a set of observable quantitaties 
which provide robust quantitative information on the network geometry. The 
lesson learned from the study of twodimensional vasculogenesis is that per- 
colative exponents are an interesting set of such observables, so we tested the 
computation of percolative exponents on simulated network structures. 

A quantitative comparison of the geometrical properties of experimental and 
computational network structures will become possible as soon as an adequate 
amount of experimental data, allowing proper statistical computation, will be- 
come available. 
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Figure 8: Cluster percolation with cell density n = 2500cells/mm'^. A: con- 
nected clusters in a realization of model (QJ. B: the largest cluster depicted in 
A percolates. 
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Figure 9: Percolation probability at varying densities A: B: the curves in A are 
collapsed according to formula ifTHl 
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In order to compute the robust statistical observables described in the paper 
one has to perform many runs of the simulation code using different random 
initial data. This, toghether with the intensive use of computational resources 
required by a three-dimensional hydrodynamic simulation on fine grids, calls for 
an efficient implementation of the computational model on parallel computers, 
as the one we presented in this paper. 

Simulations of blood vessel structures can in principle present practical impli- 
cations. Normal tissue function depends on adequate supply of oxygen through 
blood vessels. Understanding the mechanisms of formation of blood vessels has 
become a principal objective of medical research, because it would offer the pos- 
sibility of testing medical treatments in silicio. One can think that the dynami- 
cal model (QJ can be also exploited in the future to design properly vascularized 
artificial tissues by controlling the vascularization process through appropriate 
signafing patterns. 
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